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The decay of kinetic helicity is studied in numerical models of forced turbulence using either an externally imposed 
forcing function as an inhomogeneous term in the equations or, alternatively, a term linear in the velocity giving rise to 
a linear instability. The externally imposed forcing function injects energy at the largest scales, giving rise to a turbulent 
inertial range with nearly constant energy flux while for linearly forced turbulence the spectral energy is maximum near 
the dissipation wavenumber. Kinetic helicity is injected once a statistically steady state is reached, but it is found to decay 
on a turbulent time scale regardless of the nature of the forcing and the value of the Reynolds number. 
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1 Introduction 

The physical properties of astrophysical turbulence are 
often studied by solving the hydrodynamic equations in a 
periodic domain with an assumed forcing function. In par- 
ticular, isotropic homogeneous turbulence is often studied 
as a proxy of turbulence in more complicated situations, 
where specific concepts and general aspects are harder to 
isolate. An important concept is that of the forward cascade 
of kinetic energy to smaller scales. This leads to a k^^/^ 
energy spectrum, where k is the wavenumber Many flows 
of geophysical and astrophysical relevance are subject to 
rotation and stratification and can therefore attain kinetic 



direct numerical simulations ( 


Borne & Orszag 1997 


Brandenburg & Subramanian 2005a 


1, as well 



as shell model calculations JChkhetiani 19961 
IDitlevsen & Giuhani 200ralbl [Stepanov et al. 2009| l 

show an approximate k~^^^ scaling for the kinetic he- 
licity, suggesting that kinetic helicity too is subject to 
a forward cascade toward smaller length scales. On the 
other hand, kinetic helicity is conserved by the quadratic 
interactions and might therefore play an important role in 
the inviscid limit. Although this is also the case in ideal 
magnetohydrodynamics (MHD), where magnetic helicity 
is also conserved by the quadratic interactions, there is a 
significant difference. In the ideal case of small magnetic 
diffusivity, magnetic helicity can only evolve on resistive 
time scales. This can have profound effects on the saturation 
behavior of large-scale dynamos ( [Brandenburg 200 T I. 



= uxu;— VP- 



-/+iy[iV(V-tt)-VxVxM+G],(l) 



As mentioned above, in a polytropic flow, the quadratic 
interactions conserve kinetic helicity, {u) ■ u), where u) = 
V X tt is the vorticity, u the velocity, and angular brackets 
denote volume averaging over a closed or periodic domain. 
This conservation property becomes evident when writing 
the Navier-Stokes equation in the form 

du 
'dt 

where P — ^u^ + his the sum of specific turbulent pres- 
sure, and specific enthalpy, h = In p. Furthermore, 
p — const is the density, and / is a forcing function. In 
the absence of forcing, / = 0, kinetic helicity is just sub- 
ject to viscous decay, because the nonlinear term tt x a; is 
perpendicular to lj, i.e., we have 

= -Ml 



(2) 
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where q = V x a; is the curl of the vorticity. In the ideal 
case, 1/ = 0, we have {u ■ u) — const. However, in fluid 
dynamics the ideal case is hardly representative of the limit 
of large Reynolds numbers, where — >^ 0. Indeed, for a 
self-similar decay of kinetic energy, the wavenumber of the 
energy-carrying eddies, k[, decreases with time such that 
vkft « const. This implies that the rate of energy de- 
cay, v{ijj'^), is essentially independent of v and hence in- 
dependent of the Reynolds number Given that {u^) is re- 
lated to the kinetic energy, which is also independent of the 
Reynolds number, we expect that the ratio 

(3) 

which is related to the Taylor micro-scale wavenumber 
fcxay, should be proportional to Re, and therefore 

fcTay ^ Re^/^ (4) 
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However, for helical flows the rate of kinetic helicity dissi- 
pation is proportional to v{q ■ uj). Thus, if we define 



(q • ■ u) 



fr2 



(5) 



we see that kinetic helicity dissipation is related to kinetic 
energy by altogether 3 wavenumber factors. If all these 
factors scale like in Eq. (|4]i, we may expect the rate of 
kinetic helicity dissipation to diverge with decreasing v 
like 1^-1/2. This is in stark contrast to the related case of 
magneto-hydrodynamic turbulence where, following simi- 
lar reasoning, the magnetic helicity dissipation converges 
to zero like r/^/^ as the magnetic dissipation -q goes to 
zero ([Brandenburg & Subramanian 20051^; see also the ap- 
pendix of B randenburg & Kapyla (2007)| for a clear exposi- 
tion of these differences. 

A problem with the simple argument above is that in 
cases of practical relevance the forcing function / usu- 
aUy breaks kinetic helicity conservation. This is partic- 
ularly evident for the so-called linear forcing model of 



Lundgren (2003) where 



f = Au 



(6) 



is a positive multiple of the velocity vector In that case we 
have 



— (u) ■ u) = 2A{u) ■ u) — 2v{q ■ u:) 
dt 



(7) 



so that {u) ■ u) could even exhibit exponential growth. An 
aim of this paper is thus to investigate to what degree ki- 
netic helicity is conserved in forced turbulence using both 
the linear forcing model and compare it with the more tra- 
ditional stochastic forcing in a narrow wavenumber band. 
Another motivation is the fact that, by analogy, magnetic he- 
licity turned out to be of crucial importance in understand- 
ing the saturation properties of helicaUy forced dynamos in 
periodic domains (see, e.g., [Brandenburg 200 1) . We study 
the kinetic helicity evolution by monitoring the response to 
adding a large-scale helical component to the flow for both 
types of forcing. 



2 Method 

In the following we present results obtained by solving the 
compressible hydrodynamic equations with an imposed ran- 
dom forcing term and an isothermal equation of state, so 
that the pressure p is related to p via p = pc^, where Cg is 
the isothermal sound speed. We thus deviate from the con- 
ceptually simpler incompressible case. The reason for do- 
ing this is that we are mainly interested in astrophysical ap- 
plications where the gas is compressible. Furthermore, we 
expect that for small Mach numbers, compressible and in- 
compressible cases become nearly identical. However, this 
is only partially true, as the zero-Mach number limit may 
not be uniform and may not commute with the small-scale, 
long-time or zero- viscosity limits. 



In the following we use an isothermal equation of state 
for a monatomic gas for which the bulk viscosity vanishes, 
so the hydrodynamic equations for p and u take the form 

dp 
dt 



— = - V • pu, 



(8) 



du _ 

-— = —u ■ Vu 

dt 



clV\np + f + p-^V ■2pvS, 



(9) 



where Sij — \{ui,j + Uj ,i) — |<5yV • u is the traceless 
rate of strain tensor, v is the kinematic viscosity, and / is 
a forcing function that is either given by the linear forcing 
model using Eq. ^ with a positive constant A, or, alterna- 
tively, it is given by a random forcing function consisting 
of plane transversal waves with random wavevectors k such 
that |fe| lies in a band around a given forcing wavenumber 
fcf. The vector k changes randomly from one timestep to 
the next, so / is 5 correlated in time. The forcing ampli- 
tude is chosen such that the Mach number. Ma = itrms/cs, 
is about 0.1. Here, Uimg = [v?)^/"^ is the root-mean-square 
(rms) velocity. We use triply-periodic boundary conditions 
in a Cartesian domain of size . The smallest wavenumber 
that fits into the computational domain is fci = 27r/L. 

For the linear forcing model we choose for the ampli- 
tude A/cgki = 0.02, while for the random forcing function, 
/ is of the form 



f{x, t) = Re{iV/fe(t) exp[ifc(t) • x + i0(i)]}, 
X is the position vector and 



/fc = (fcxe)/v'fc2-(fc-e)2, 



(10) 



(11) 



where e is an arbitrary unit vector not aligned with k; note 
that l/fcp = 1. The wave vector k{t) and the random phase 
— TT < (t){t) < TT change at every time step. For the time- 
integrated forcing function to be independent of the length 
of the time step 5t, the normalization factor N has to be 
proportional to St^^^^. On dimensional grounds it is chosen 
to be = foCs{kfCs/6ty/^, where /o is a nondimensional 
forcing amplitude, which is chosen to be /o = 0.02. For 
the monochromatic ally forced simulations we choose fcf = 
1.5 fci. For the linear forcing module, on the other hand, we 
compute fcf is the integral scale from the resulting kinetic 
energy spectrum (see below). 

Our simulations are characterized by the value of the 
Reynolds number. 



Re = Urrns/vkf 



(12) 



where the wavenumber fcf is either the forcing wavenum- 
ber in the case of monochromatic random forcing, or it is 
evaluated as the wavenumber of the energy-carrying eddies. 



fc; 



J k-^E{k)dk I J E{k)dk, (13) 

with E{k) being the kinetic energy spectrum, which is here 
normalized such that 

/•OO 



E{k) dfc 



(14) 
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Fig. 1 Kinetic energy spectra compensated with 
fc5/3g-2/3 Jqj. ^ range of different Reynolds numbers and 
numerical resolutions up to 512^ mesh points. 



It turns out the in the former case of monochromatic ran- 
dom forcing, the definition (flJl i agrees well with the a pri- 
ori chosen forcing wavenumber We plot energy spectra as 
a function of k in units of the dissipation wavenumber, 

fed = (15) 

Here, e = i^{<^'^) is the rate of energy dissipation. We also 
analyze spectra of kinetic helicity, F{k), which are normal- 
ized such that 



POO 

/ F{k)dk = {lo -u) 
Jo 



(16) 



We note that the kinetic helicity spectrum can have either 
sign and its modulus obeys the well-known realizability 
condition, \F{k)\ < 2kE{k) (IMoffatt 19691 1. 

In the statistically steady state, the kinetic helicity fluc- 
tuates around zero. We consider the evolution of kinetic he- 
licity after adding an instantaneous finite perturbation in the 
form of a Beltrami field. 



It — 5- tt + ecs(sin fcz, cos fcz, 0), 



(17) 



where e is a nondimensional input parameter. In all cases we 
choose k = ki, i.e., we perturb the system with a wave on 
the scale of the domain. We express time in units of turnover 
times, so that tu^a^Jtf is a nondimensional time. 

The simulations have been carried out using the PENCIL 
Code which is a high-order finite-difference code (sixth 
order in space and third order in time) for solving the com- 
pressible hydrodynamic equations. 

3 Results 

3.1 Linear forcing model 

In Fig. [U we show compensated spectra of kinetic energy 
for different values of the Reynolds number. The spec- 
tra collapse onto each other in the dissipation range if the 



wavenumber is scaled with k^- With increasing Reynolds 
number the spectra begin to sketch out a continuation to- 
ward smaller wavenumber k. Note that in the linear forc- 
ing model, energy is being injected at all length scales 
and not just at the largest scale of the domain. This is 
a property that may also be responsible for the fact that 
the effective driving scale tends to be smaller than in oth- 
erwise equivalent monochromatic ally forced simulations 
dRosales & Meneveau 2005l l. 

In the statistically steady state, the kinetic helicity fluc- 
tuates around zero. However, after adding a finite perturba- 
tion in the form of Eq. ([17), the helicity begins to grow for 
about 5 turnover times, (5iMimsfcf ~ 5. After that time there 
is a systematic decline of the kinetic heUcity. By inspecting 
the results for different values of Re, we see that the time 
it takes for the kinetic helicity to relax to previous levels 
becomes longer as the Reynolds number is increased from 
44 to 125; see Fig. |2] However, for Re = 300, which is 
our largest Reynolds number for which we have performed 
experiments with added Beltrami fields, the decline of he- 
licity has happened on a similar time scale as for Re = 70. 
This suggests that there may not be a systematic Reynolds 
number dependence of kinetic helicity decay. 

It should be noted that during the first 2-3 eddy turnover 
times after adding the Beltrami field perturbation, the ki- 
netic helicity shows an exponential increase; see also Fig. [3] 
This is connected with the fact that in the absence of any 
other effective damping, Eq. (|7]i would imply an exponen- 
tial growth of {ijj ■ u) (X. exp Xt with growth rate 



\ = 2{A- vkl^), 



(18) 



http://pencil-code.googlecode.com/ 



where k^fi quantifies the involvement of high wavenum- 
bers in the expression for the kinetic helicity dissipation; 
see Eq. (|5]l- 

One would expect strong perturbations to survive this 
exponential growth for longer, but this is not the case, as 
is demonstrated in Fig. [3] where we show the evolution of 
the kinetic helicity for different values of e after adding the 
Beltrami field perturbation. The reason for the subsequent 
decay of magnetic helicity lies in the fact that fcoff scales 
like fceff oc Re^/^, so the rate of kinetic helicity dissipation 
remains significant even in the limit Re oo, i.e. — ^ 0. 
This is quantified in terms of k^tf, whose scaling will with 
Reynolds number will be considered below. 

3.2 Monocliromatic random forcing 

In order to see whether the results presented above are a 
special consequence of the linear forcing model, we now 
perform simulations using the more traditional monochro- 
matic forcing in a narrow wavenumber interval. In Fig. |4] 
we show energy spectra for different Reynolds numbers. 
The results suggest that the decline in spectral power toward 
the smallest wavenumbers seen in Fig.[T]for the linear forc- 
ing model is now absent. In other words, while in Fig.|4]we 
clearly see that the compensated energy spectrum is flat, this 



www.an-joumal.org 



©2012 WILEY- VCH Verlag GmbH & Co. KGaA, Weinhdm 



198 



A. Brandenburg & A. Petrosyan: Kinetic helicity decay 



Re=44 




-40 -20 20 40 60 80 100 



Re = 70 



A 
V 

■sT 
\ 

A 

3 

V 




-40 -20 20 40 60 80 100 



Re=125 




-40 -20 20 40 60 80 100 



Re=300 




-40 -20 20 40 60 80 100 
(i-t.) WrniA 



Fig. 2 Evolution of the normalized kinetic helicity for 
different values of Re after adding a Beltrami field pertur- 
bation at t = t, . 
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Fig. 3 Evolution of the kinetic helicity for different val- 
ues of e after adding the Beltrami field perturbation att — 
U. 
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Fig. 4 Kinetic energy spectra compensated with 
/i;5/3g-2/3 fQj. ^ range of different Reynolds numbers and 
numerical resolutions using monochromatic forcing. 

is not the case in Fig.[T]for the linear forcing model. How- 
ever, there is still an uprise near 0.1 fed that one may gen- 
erally associate with the bottleneck effect dFalkovich 19941 
IDobler et al.^OOST l. 

Next, we study the effect of adding a helicity perturba- 
tion also in this case. Figure|5]gives time series for three val- 
ues of Re. There is no evidence for a prolonged relaxation to 
zero. The reason for this could be that a helical wave cannot 
interact with itself; see a corresponding discussion follow- 
ing Eq. (11) of [Kraichnan ( 1 973)| This also suggests that 
kinetic helicity conservation in the inviscid case, v — 0, is 
of no relevance to the inviscid limit, v 0, in which case 
the kinetic helicity dissipation diverges. This behavior was 
less clear in the previous case with the linear forcing model. 
This might either be a matter of coincidence, but it could 
also be a consequence of the linear forcing model which 
has exponentially growing solutions. 

The time series in Fig. |5]reveals another interesting as- 
pect in comparison to Fig. |2]for the linear forcing model in 
that the level of fluctuations of (u) ■ u) is generally larger for 
the monochromatic forcing function than for the linear forc- 
ing model. Furthermore, the time series show much stronger 
short-term fluctuations while for the linear forcing model 
the time traces of {lj ■ u) are smoother. 

In Table [T] we summarize integral and dissipation 
wavenumbers as well as the normalized energy fluxes for 
both linear and monochromatic random forcings. In wind 
tunnel turbulence one usually expresses the energy flux in 
units of a quantity = u'^ / L, where u' = Mrms/\/3 
is the one-dimensional rms velocity and L — Sir/Aki is 
the customary definition of the integral scale. The stan- 
dard result of Ce ~ 0.5 dPearson et al. 20041) corresponds 
then to e ~ 0.04:k{u^^g. Comparing the normalized en- 
ergy fluxes for linear and monochromatic random forcings 
we see hardly any differences. This suggests that the nature 
of the forcing in hydrodynamic turbulence might not be of 
great qualitative importance, although it is still possible that 
the botdeneck effect dFalkovich 19941 IDobler et al. 20031 ) 
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Table 1 Summary of the normalized characteristic wavenumbers k{ — kf/ki, kd = kd/ki, /cTay — fcTay/fci, and 
fceff = kcs /ki, for linear and monochromatic forcings. The numerical resolution is given in the second column. 
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near the dissipative scale might be connected with the na- 
ture of the forcing at large scales (IDavidson 2004l l. 

In Table □ we also list the values of fcxay and fceff , as ^ig. 6 (onUne colour at: www.an-journal.org) VisuaUza- 

defined in Eqs © and ©. Both wavenumbers are clearly ^lon of the logarithmic density (proportional to logarithmic 

proportional to Re^/'\ as can be seen from Fig.|71 where we Pressure fluctuations) on the periphery of the domain for hn- 

plot the Reynolds number dependence of the ratios fcxay /fci ^"'^ forcings. In both cases we have Re « 600 

and fceff/^i for linear and monochromatic forcings. ^ resolution of 512 mesh points. 

To compare the linearly and monochromatic ally forced 
turbulence simulations in real space, we present in Fig. |6] 

visualizations of the logarithmic density on the periphery convenient scalar quantity characterizing the pressure fluc- 

of the computational domain. The logarithmic density is a tuations resulting from the Reynolds stress. These visualiza- 
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Fig. 7 (online colour at: www.an-journal.org) Scalings 
of fcxay/fci (solid Hnes) and kcs/ki (dashed lines) with Re 
for linear (filled symbols) and monochromatic (open sym- 
bols and red lines) forcings. The dotted line indicates Re ' 
scaling. 
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Fig. 8 Comparison of compensated kinetic energy and 
helicity spectra for the linear forcing model (solid line) 
and the monochromatic random forcing function (dashed). 
In those patches where the kinetic helicity is negative, the 
modulus of the value is plotted as a dotted line (red for the 
linear forcing model and blue for the monochromatic ran- 
dom forcing function). 



tions look rather similar in the two cases, suggesting that the 
slight difference seen in the power spectra in Figs.[T]and|4]is 
hard to discern in such images. In fact, upon closer compar- 
ison of just the cases with the highest resolution we see that 
both kinetic energy and kinetic helicity spectra agree almost 
perfectly (Fig. |8]l, except that in the linear forcing model 
the kinetic energy spectrum has dropped significantly at the 



smallest wavenumber, which is not the case for monochro- 
matic random forcing. 

4 Conclusions 

The linear forcing model is based on the driving of turbu- 
lence by a linear instability instead of a forcing function that 
is independent of the flow. This type of forcing might be 
more physical, because it does not change abruptly and de- 
pends on the local flow properties. Nevertheless, both types 
of forcing are Galilean invariant. This would change if the 
flow-independent monochromatic forcing were no longer S 
correlated in time. A disadvantage of linearly forced tur- 
bulence is that energy injection occurs at all wavenum- 
bers. Our simulations show that the energy spectrum is 
perhaps slightly shallower than with the flow-independent 
monochromatic forcing function. Part of this is explained by 
flie bottleneck effect dFalkovich 19941 IDobler et al. 20031 1. 
and that the energy spectra compensated by k^^^ show a 
stronger uprise toward the dissipation wavenumber. It is un- 
clear whether this result would persist at larger resolution. 

The linear forcing model has the interesting property of 
amplifying not only kinetic energy, but also kinetic helic- 
ity. Indeed, at early times, just after having injected kinetic 
helicity into the system, the coherent helical part continues 
to increase exponentially, but the flow soon breaks up into 
smaller eddies, giving rise to enhanced effective dissipation. 
For practical applications, it should be noted that the lin- 
ear forcing model has the disadvantage that the initial expo- 
nential growth of kinetic energy and kinetic helicity might 
prevail for too long. This will be the case when the initial 
random perturbations are too weak to perturb the flow suf- 
ficiently. In that case, the kinetic energy would quickly in- 
crease to large values without producing three-dimensional 
turbulence. 

With regards to geophysical and astrophysical applica- 
tions we can say that in a turbulent system, kinetic he- 
licity is no longer a conserved quantity, even though it 
would be if = were strictly true. The latter require- 
ment is of course not really possible in a turbulent sys- 
tem, because kinetic energy would then accumulate at the 
smallest possible scale resolved within the hydrodynamics 
framework and kinetic energy would not be able to decay, 
which is unphysical. While this should not be surprising, 
it is important to remember that this is quite different in 
the case of magnetohydrodynamics, where magnetic helic- 
ity dissipation really does go to zero in a turbulent system 
- even for finite (but small) values of the magnetic diffu- 
sivity. At the same time, magnetic energy dissipation does 
stay finite and is able to accomplish magnetic reconnec- 
tion on the smallest resolved scales of the turbulent cascade 



( |Galsgaard & Nordlund 1996} ILazarian & Vishniac 1999l l. 

This paper has also shown that, regardless of the 
nature of the forcing, there are fairly strong helicity 
fluctuations. They appear to be coherent over many tur- 
bulent eddy timescales. One may wonder how generic 
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such fluctuations are and if such fluctuations could be 
relevant for say the incoherent dynamo effects that have 
been investigated by several authors in recent years 
dVishniac & Brandenburg 1997' Sur & Subramanian 20091 
iHeinemann et al. 201 ll Mitra & Brandenburg 2012[ 



IRichardson & Proctor 201(jb . It will therefore be interesting 
the associate the kinetic helicity fluctuations with those of 
a, which have already been determined in simulations of 
turbulent shear flows ( Brandenburg et al. 2008[ l. 
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